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Abstract 

> 

■ Starting with an exact and simple geodesic, we generate approximate geodesies by summing up 
higher-order geodesic deviations within a General Relativistic setting, without using Newtonian 

■ and post-Newtonian approximations. 

We apply this method to the problem of closed orbital motion of test particles in the Kerr 
metric space-time. With a simple circular orbit in the equatorial plane taken as the initial geodesic 
we obtain finite eccentricity orbits in the form of Taylor series with the eccentricity playing the 
role of small parameter. 

The explicit expressions of these higher-order geodesic deviations are derived using successive 
systems of linear equations with constant coefficients, whose solutions are of harmonic oscillator 
type. This scheme gives best results when applied to the orbits with low eccentricities, but with 
arbitrary values of (GMJ Rc 2 ), smaller than 1/6 in the Schwarzschild limit. 



PACS number (s): 04.25.-g, 04.25.Nx 



1 Introduction 

In two recently published articles, Jfl] and PJ, the approximate low-eccentricity relativistic trajectories 
of planets with small mass m (as compared to the central body's mass M) have been constructed 
in Schwarzschild and Reissner-Nordstr0m space-time metrics. In the latter case, the motion of 
electrically charged particles have been investigated, too. 

The two-body problem in General Relativity has been the object of many excellent studies; one 
of the first checks of this theory has been the very precise value of the perihelion advance calculated 
by Einstein || for the planet Mercury. The calculus was based on the solution of the geodesic 
equation in Schwarzschild's metric, using the first integrals; the solution was obtained in the form 

*e-mail: coliste@ccr.jussieu.fr 
^e-mail: leygnacOlapp . in2p3 . f r 
'"e-mail: rk@ccr.jussieu.fr 



1 



of a quadrature, with the proper time r expressed as a quasi-elliptic integral. Such an integral 
can not be evaluated analytically; instead, Einstein has developed the integrand into a power series 
with respect to the small parameter GM/rc 2 , which led to simple integrations that could be easily 
performed. The approximate formula for the perihelion advance after one revolution is then 

6nGM 6ttGM . 2 4 6 
aJT^e 2 ) = q (1 + e + e + e + ...) (1) 

where a is the major semi-axis and e the eccentricity of the orbit, and G stands for Newton's gravi- 
tational constant divided by c 2 , so that in this notation the quantity GM/r becomes dimensionless. 
We shall adopt this notation from now on, equivalent of using the units in which c = 1. The above 
approximation is acceptable only if the value of the parameter GM/a is negligibly small. 

In our previous articles H and |lj we have proposed an alternative way of determining the value of 
perihelion advance and finding an explicit (parametrized by the proper time) form of trajectory and 
the law of motion as a series of successive approximations, without supposing that the ratio GM/a (or 
equivalently, the ratio v 2 /c 2 ) is very small, or using the Newtonian approximation. Instead, we start 
from a very simple particular solution of geodesic equation in Schwarzschild (or Reissner-Nordstr0m) 
metric: a perfect circular orbit along which the small mass m is advancing with a constant angular 
velocity. It is very easy to check that such a motion is a geodesic curve in the aforementioned 
space-times. The fact that all geometrical quantities, such as the Christoffel coefficients, or the 
components of Riemann's tensor, take on constant values on this trajectory, leads to a particularly 
simple form of the geodesic deviation equations: they reduce themselves to a system of second- 
order linear differential equations with constant coefficients, and the solution is just a collection of 
harmonic oscillators. 

Here we apply the same method to the case of the axially-symmetric Kerr metric. Although the 
motion of test particles along closed orbits in Kerr's metric has been analyzed in a very exhaustive 
manner in many papers our method gives results in an explicit form, and is very well adapted 

for computer-based calculations. In the case of orbits with low eccentricity it converges very quickly 
even for the non-negligible values of the ratio GM/r. 

In the following Sections we shall briefly recall the essential features of our approximation method 
best suited for higher-order deviations. Using the Kerr space-time, we choose a circular orbit in the 
equatorial plane as the initial geodesic. Then, the first, second and third deviations are obtained, the 



latter one with the help of the Poincare's method [11| enabling us to obtain higher-order corrections 



to the basic frequencies. The explicit form of the perihelion advance in the field of Kerr metric 
displays interesting features as a combined result of the influence of two essential parameters, the 
mass M and the angular momentum density a of the central body. 

In the last section, we discuss the physical content of the results and consider some future 
applications of higher-order deviations, including the effects of finite mass m and internal spin 
(angular momentum) of the planet. In contrast with our previous article we shall not consider 
here the problem of gravitational radiation; it will be left to a detailed future work. 



2 Geodesic deviations using small deformations 

The previous article [JIJ was based on the deviation vectors n^, b^, hP and their deviation equa- 
tions. However, as the order of the deviation increases, it becomes harder to calculate the deviation 
equations for the deviation vectors. 

So, for our purpose here, which is the effective calculation of deformations of circular orbits 
in Kerr metric, we need the explicit coordinate-dependent expressions for the deviations that we 
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shall add to given functions of proper time s which define the relativistic trajectory and the law of 
motion. This is why we need to consider an alternative approach, which deals with small deviations 
of arbitrary order of coordinate functions, thus deforming the trajectories directly. 
Consider an infinitesimal deformation of the geodesic curve x^ (s): 

xM( a ) ^^(s) = x fl (s) + 5x> 1 (s). (2) 

Suppose that we want the new curve x^(s) to satisfy the geodesic equation, too: 

d x^ ii , _ . , , dx dx _ , _ >, 

— + rUx u ) — — = 0. 3 

ds 2 xp ds ds 
Expanding the Christoffel coefficients into power series of Sx^, 

= + to* r$,0O + 1 fe* dx- d 2 aT I*,(s") + ... (4) 

Substituting the expressions (0) and (0j) into Eq. (|3|) and expanding in consecutive powers of devia- 
tions Sx^, we get in the zeroth-order the initial geodesic equation @ satisfied by x p (s); collecting 
then all the terms linear in 5x p and their derivatives, we get 

d ^f + 2 + (d, rj,) « Vfa* = (5) 

which coincides with the first-order deviation equation of the Ref. |jj if we replace the vector n p 
by the infinitesimal deviation 5x p . The geometrical meaning of this equation is now very clear: 
it gives the conditions to be satisfied by infinitesimal functions en p (s) = 5x p (with e being an 
infinitesimal parameter) defined along a given geodesic curve x x (s), in order to ensure that the new 
curve, infinitesimally close to it and defined by x x (s) = x x (s) + en x (s) , is also a geodesic one, up 
to the first order in e. 

The fact that n p is a vector is in agreement with the transformation properties of infinitesimal 
deviations 6x^, which under an arbitrary change of coordinates x p = x p (y p ) transform as 

^ = |aV (6) 

up to higher-order terms, neglected at the linear approximation level. 

However, the higher-order terms in the expansion are still there: collecting all the second-order 
terms in Sx^ from the Taylor expansion of Eq. @, we get 

and there is no reason for it to vanish even if the deviations 5x ll (s) satisfy Eq. (||). The vanishing 
of expression (|jj) would impose too many conditions, a priori incompatible with Eq. (||) on the 
same set of functions 5x p (s), and we need extra degrees of freedom if we want to cancel also all the 
second-order terms. 

This means that from the very beginning, infinitesimal deviations of higher-order must be intro- 
duced: 

x p [s] = x p (s) + 5x p (s) + j tfV (s) + y <5V (s) + . . . (8) 
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so that two new second-order terms will add up to the expression ml), namely 

^P + ^^,r; u V ( 9) 

and the sum of all these terms represents a new set of second-order differential equations imposed 
on the independent functions 5 2 x p , which may be solved after we insert the solutions for it^(s) and 
Sx^ts) obtained previously. 

At this point, another problem arises: not only the coefficients in these differential equations 
are not covariant objects, but also the functions S 2 x >M do not behave as vectors under the change of 
coordinates. As a matter of fact, they will mix up with terms quadratic in 5x^ as follows: 

*V = 6(5xn = 5 2 y x + 5y X 6y» . (10) 

This non-homogeneous transformation law suggests that we can introduce a covariant quantity D 2 x^ 
defined as 

D V = $V + It p 5x x 5x p . (11) 

Defining the infinitesimal vector b p as e 2 = D 2 x p and expressing the Taylor expansion (|8|) in terms 
of n M and b^ as 

+ \f 2 (> - n x n p ) +... (12) 



2! V" "V' 

and requiring the geodesic equation for x^ 1 to be satisfied up to the second order in e, we arrive at the 
same second-order deviation equation of the Ref. Q satisfied by b^ that is non-manifestly covariant 
and equivalent to the manifestly covariant equation. 

In Ref. Q, the non-manifestly covariant deviation equations could be easier obtained using 
this Taylor expansion approach than transforming the manifestly covariant deviation equations. In 
practical calculations the non-manifestly covariant equations turn out to be of much use, i.e., it is 
easier to obtain the solutions of and b^ using them instead of the manifestly covariant equations. 

Similar corrections are needed to define the higher-order deviations, like 

e 3 h p = D 3 x^ = 5 3 x p + 3T p a 5x u 5 2 x a + (d x T p a + T p p r£ a ) 5x x 5x u 5x a , (13) 
so the real coordinate deviation 5 3 x >M reads 

^ = £ 3^ _ n p b a _ ^ T ^ _ ^ T P^ n A n , n1 . (14) 

A study of higher-order differentials and their covariant generalizations can be found in recent papers 

But there is even another easier form to calculate higher-order geodesic deviations : we indeed 
need 5 n x p to obtain the geodesic x p , and the differential equations for 5 n x p are simpler than their 
counterparts n^, 6 M , W, etc. For example, requiring again the geodesic equation for x M to be satisfied 
up to the second order in e, the following second-order deviation equation for 5 2 x p is obtained 

+ (d p Fl) «V«V + 2T^u xd -^f = 

~ 2 r v^^r " 4(a - Txp) uX 5x °^f - {d » d ° Txp) uX uP 5x ° 5x "> (15) 
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where we see that the l.h.s. is unchanged, but the r.h.s. has only 3 terms instead of 10 found in the 
non- manifestly covariant second-order deviation equation for (see Ref. [||). 

The non-manifestly covariant third-order deviation equation for hP is not shown here, but has 
60 terms in the r.h.s., while the third-order deviation equation for S 3 x fJ ' has only 7 terms: 



ds 2 y 1 Aa ' ds 

„ d5 2 X X ddxP ( a d5x X d5 X P A 2 A r CT *>V\ 

_6r L — ; ; oOffrL [ox — ; hit ox — \- u dx — ; — 

Xp ds ds v Xp> \ ds ds ds ds J 

-3(0„ 8 a F p p ) u x 5x u (2 Sx a ^- + u" 5 2 x< 7 ^j -{d T 3 V d a T p p ) u x u" 5x a 8x u 5x T . (16) 

The fourth-order deviation equation for 5^x p has 15 terms, and the fifth-order deviation equation 
for S 5 x p has 26 terms. We have developed a symbolic computer program to calculate n^-order 
deviation equations for S n x p . 

The non-manifestly covariant geodesic deviation equations are well suited to deriving successive 
approximations for geodesies close to an initial one. Starting from a given geodesic x^ (s) we can 
solve Eq. (|5|) and find the first-order deviation vector 5x p (s). Now, with (s) and Sx^ (s), the 
system ( |i~5| ) can be solved and we obtain the second-order deviation 5 2 x p ' (s).Then, using u M (s), 
5x^ (s) and 5 2 x Pj (s) into the system (|i~6| ) the third-order deviation 5 3 x p (s) is calculated, and so 
forth. 

The literature about geodesic deviations includes a rigorous mathematical study of geodesic de- 
viations up to the second-order, as well as geometric interpretation, but using different derivation, 
presented in |]14|. Also, a Hamilton-Jacobi formalism had been derived in |]15||, which was applied 



to the problem of free falling particles in the Schwarzschild space-time [16]. Anyway, the resulting 
expressions are not well optimized for successive calculations of higher-order geodesic deviations. 
Interesting effects resulting from the analysis of first-order geodesic deviations of test particles sus- 
pended in hollow spherical satellites have been discussed in |fT7f| . 

3 Circular orbits in the Kerr metric 

We choose as initial geodesic a circular orbit in the axisymmetric gravitational field created by a 
massive body with rotation, i.e., in the Kerr metric. This metric and their circular orbits have been 
studied in several papers and books ||-|10|. 



The gravitational field is described by the line-element (in natural coordinates with c = 1 and 
G=l) 



o p 2 n 00 2Mr — p 2 AMra Q sin 2 (/ , .-, , .-, , ., 

ds 2 = ^-dr 2 + p 2 d9 2 + TT^dt 2 7 ^sin 2 9dtd(t)+ — 5— -Aa sin 6 + (r +a ) ) d<j) 17 



with A = r 2 + a 2 — 2M r and p 2 = r 2 + a 2 sin 2 #, where M and a = 4j are the mass and the angular 
momentum density of the central body rotating in opposite 4> direction. 

The circular orbit of radius R in the equatorial plane (which is a geodesic in the background 
Kerr metric) is described by a simple 4-velocity vector: 

dr dO 

u r = — = 0, u 9 = — = 0, 
as as 
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V R R3/2 

because r = R = const., 9 = tt/2 = const., so that sin 9 = 1 and cos 9 = 0. The angular 
frequency of this circular motion is u c /u , which decreases if the central body increases the rotation 
in the opposite direction. 



4 First-order geodesic deviation around a circular orbit 

Likewise the calculations for the Schwarzschild case in Ref. Q, we use the non-manifestly first-order 
deviation equation for 5x^, Eq. (||). The vector u M of a circular orbit, calculated before, is used, 
providing constant components that yield simple equations for the components Sr, 69, 5<p and 5t (or 
n r , n 6 , r& and n t ), as we can see in a matrix form: 



to 2 i m 22 m 2 3 m 2 4 
?7i3i m 32 m 33 m 34 



\ 


/ Sr \ 
59 








<l\ 




5<p 







/ 


\St) 




W 



(19) 



where the matrix elements are: 
3M/ 2 



ds 2 



R 3 fx 



77112 = 0, 77113 



1 



m 3 i 



m 4 i 



m 2 i = m 23 = m 24 
2^ / 4 d 

2M / 5 d 
" ^ 2 /avTTds' 



#3/2 
d 2 



using the functions: 



°' m22 = d? + i^X' 



mi4 



"132 = "T-34 = 0, m 33 



m 42 = m 43 = 0, m 44 



~~ ds 2 ' 

£_ 

ds 2 ' 



2M f 2 d 
R 2 vTTds 



, (20) 

(21) 
(22) 
(23) 



3M\ 



+ 



2ax/M 



R 3/2 
1 



2M\ 



+ 



R 2 ' 



4aVM 3a z 

#3/2 + #2 



/4 = 1 - + 



2M 
~R 



aVM 



#3/2 



2aVM a z 
~W 2 ~ + R 2 ' 



The harmonic oscillator equation for n 6 



io 



M h 



59 has an angular frequency ojq: 



R?l 2 V /i R?l 2 



1 



4aVM 



3a 2 
~W 



+ 



2aVM 
#3/2 



(24) 



(25) 
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One possible choice of solution is 



n 



69 



-n cos(uigs). 



(26) 



In the Schwarzschild limit (a — ► 0), ojq = lj c , so in this case we can neglect this solution (n® = 0) 
because the new plane of orbit is a new one inclined, or just a change of coordinate system. 



Using the differential equation for n r 
the harmonic oscillator equation 



with the characteristic frequency 



8r we can eliminate the derivatives of 5<fi and St, yielding 



d 2 5r 
ds 2 



+ u 2 5r = 0, 



(> 



6M\ 
R ) 



+ 



8aVM 

R 3/2 



3a 2 
~W 



M 

2 




where: 



6M\ 



+ 



8aVM 



3a 2 



We shall choose the initial phase to have (with n r > 0): 

n r = 5r = —tiq cos(o;s) 

so the perihelion occurs when s = 0. 

The calculation of 5(f) and 5t is now simple: 



n 



5(f) = Uq sin(u;s), 
6t = n\ sin(ws), 



where the amplitudes depend on n r : 



nk = 



no = 2no 



2n r f A 




(27) 

(28) 

(29) 
(30) 

(31) 
(32) 

(33) 
(34) 



Adding this first-order deviation to the circular orbit, the new trajectory and the law of motion 
are given by 



R — n r cos(ws), 



7T 



t 



- - n cos(ujgs), 
uo c s + Uq s'm(uj s 



rin sinfw s) 



(35) 
(36) 

(37) 
(38) 



and this solution is a geodesic up to the first-order in e. It is important to note once again that the 
coefficient n r , which also fixes the values of the two remaining amplitudes, n-g and n^, defines the 

n r 

size of the actual deviation, so that the ratio becomes the dimensionless infinitesimal parameter e 
controlling the approximation series with consecutive terms proportional to the consecutive powers 

of § (o4). 
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One easily checks that the Schwarzschild limit (a — ► 0) of the solution above yields the results 
of the Ref. including the perihelion advance and the generalized epicycle [18|, where we identify 
the major semi- axis a with R and the eccentricity e with : 



r(i) 



o(l 



1 + e cos(u>o t) 



a [1 — e cos(u>o i)] 



(39) 



But it is more interesting to show the perihelion advance for the Kerr case, even using the first-order 
deviation, 



2vr 



1 



540 ttM 5 / 2 

+ RV 2 + 



'6vrM 27ttM 2 135 ttM 3 \ 
~ R~ + R 2 + ft 3 + -J 

\ , ( 3 vr 75 vrM 1845 ttM 2 



'871^ + 72vrM 3 / 2 



i? 2 



+ 



i? 3 



+ 



2i? 4 



+ 



^3/2 
+ ... 



^5/2 



(40) 



because the Kerr parameter a appears with positive and negative coefficients, i.e., the angular 
momentum density o can increase or decrease the perihelion advance. Note that the post-Newtonian 
limit matches the Eq. (||) for small eccentricities. 

Despite the limitations of the first-order geodesic deviation, we have already obtained a genera- 
lized perihelion advance valid for high values of ^ and a, but low values of the "eccentricity" e (or 
The high-order deviations will, for example, allow the calculation of Aip for higher values of 



R 



the "eccentricity" 



^0 

R ■ 



5 The second-order geodesic deviation 



Inserting the complete solution for the first-order deviation 5x^ = n^, Eqs. (fffi), (j3C|)-(|3"2|) into the 
second-order deviation equation for 8 2 x fl (|i~5|), we find the same matrix with a new non- homogeneous 
vector C M : 

/mn mi2 mi 3 m u \ 
m 2 i m 22 m 2 3 m 24 



m 3 i 
m 4 i 



m 32 
m 42 



"133 

m 43 



m 34 
m 44 



/5 2 r\ 
' 5 2 6 * 

<5 2 </> 



(41) 



where we have put into evidence the common factor e 2 , which shows the explicit quadratic depen- 
dence of the second-order deviation 5 2 x fl on the first-order deviation amplitude n r (or n®). The 
functions C r ,C e , and C* are expressions depending on M, i?, a, and on the functions sin(2u;s), 
cos(2u;s), sin(2u;6»s), cos(2ujgs), cos [(uj — ujg)s] and cos [(uj + wg)s]: 



C r = Cq + C 2r cos(2cjs) + C r 2e cos(2u e s), 

C e = Ci cos[(cj - uj e )s] + C e + cos[(u + w e )a], 
= Ct sin(2u;s) + C| e sin(2u e s), 
C* = C 2r sin(2u;s) + C\ e sm(2uj e s). 



(42) 

(43) 
(44) 
(45) 

The solution of the above matrix for 5 2 x >1 (s) has the same characteristic equations of the matrix 
( |T9| ) for &r M (s) = n^(s), and the general solution containing oscillating terms with angular frequency 
uj and ujg is of no interest because it is already accounted for by n^(s). But the particular solution 
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includes the terms linear in the proper time s, constant ones, and the terms oscillating with angular 
frequency 2u, (u> — ujq) and (u + u>g): 



5 2 r = 5 2 tq + 5 2 T2 cos(2u;s), 
5 2 = 5 2 6_ cos[(w - u B )s] + 5 2 6 + cos[(w + u )s], 
8 2 (f> = (5 2 (j)o)s + 5 2 (f>2r sin(2a;s) + 5 2 4>28 sm[2ujgs], 
5 2 t = (S 2 to)s + 5 2 t2r sin(2ws) + 5 2 t2o sin [2^0 s]. 



(46) 
(47) 
(48) 
(49) 



The constants 5 2 ro, 5 2 4>o and 5 2 to depend on two arbitrary constants, so we can choose the initial 
conditions of the differential solutions so that the constants <5 2 ro and <5 2 </>o are null, and <5 2 io is 
simplified. The Appendix 1 shows the explicit values of the above coefficients. 

In the Schwarzschild limit, the solution for the second-order geodesic deviation 5 2 x fl (s) is: 



5 2 r 



K) 2 ( 



1 



7M \ 
R ) 



5 2 t 



R (l _6M 

K) 2 



cos(2ws), 5 2 9 = 0, 5 2 



_ 2 K 



R 2 







6M\ 3 / 2 
R 



R 



1 + 



M 



R(i-^r 



1 _ 3M 
1 R 



-s + 




15M , UM 2 
R T 



2M 
R 



) 2 



6AfA 3 / 2 
R J 



sin(2a;s) 



(50) 



(51) 



The second-order deviation S 2 x^ computed in the Ref. JxJ used another choice for the constants 
5 2 ro, 5 2 4>q and 5 2 to, i.e., equivalent to different initial conditions for the initial geodesic. So, using 
the initial conditions chosen above, the comparison with an ellipse in the Schwarzschild limit also 
gives more compact results. 

The trajectory described by x^ including second-order deviations is not an ellipse due to the 
General Relativity effects of but we can match the perihelion and aphelion distances of the 
Keplerian, i.e., elliptical orbit, with the same perihelion and aphelion distances of the orbit described 
by x M : 



R 



) 



2R 



(l-Sjf 



2*5 (l - f 



2i?( 



1 



6M\ 
R J 



K) 



i 



7M \ 
R ) 



R 



r\3 



+ o 



K) 

R 3 



(52) 



(53) 



The shape of the orbit described by r(<p) can be obtained from ip(s), then s(ip) by means of successive 
approximations beginning with los = -fj-tp, and s is replaced in r(s) giving r(ip) up to the second 

order in 



??-, 



- = 1 - -H cos _ (p + 



R 



R 



n- 



R 



-1 + 



V COS (fi 

1 _ 6M\ 



+ ... 



(54) 



The exact equation of an ellipse is obtained in the limit ^ 

R 



0, up to the second order in e = Uq/R: 



(l - f e 2 



1 + e cos (p 1 + e cos ip 



R 



1 — e cos ip + e 



1+2 cos2 </? ) + ... 



(55) 
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In the ellipse equation ( |39| ) we have ro = a(l — e 

r(i 



so 



3 2 
2 e 




(56) 



These values for a and e agree with Eqs. (|52|)-(53). 

In order to improve the comparison of the perihelion advance in the post-Newtonian limit with 
Eq. (H), the Aip should include ^ terms, which is not yet the case using second-order deviations 
due to the imposed initial conditions. 



6 Third-order deviation and Poincare's method 



Using the solutions for 8x^ = and 5 2 x^ into third-order deviation equation for 5 3 x >1 (16), we 
again find the same matrix with a new non- homogeneous vector L>^: 



/mn rn 12 m 13 m M 
m 2 i m 2 2 rn 2 3 m 2 4 
m 3 i m 32 m 33 m 34 



\ 


( 5 3 r \ 




(° r \ 




1 5 3 e 1 














) 


\6H) 




KD*) 



(57) 



where the common factor e 3 shows the explicit cubic dependence of the third-order deviation 5 3 x 11 
on the first-order deviation amplitude rig (or rig). The functions D r ,D s , and D l are expressions 
depending on M, R, a, and sin and cos functions of los, loqs, 3uj, 3u>qs, (lo — 2loq), (uj + 2ujq), (2lo — loq) 
and (2u> + loq): 

D r = D\cos{los) + D 3 cos(3ws) + D r _ cos[(w - 2uj e )s] + L>+cos[(w + 2u e )s], (58) 

D e = D[ cos(uj e s) + D e 3 cos(3uj e s) + D e _ cos[(2w - u e )s] + D e + cos[(2w + u e )s], (59) 

= Df cos(ujs) + Z?3 cos(3ujs) + Dt cos[(w - 2uj 9 )s] + L>J cos[(w + 2uj 9 )s], (60) 

D* = D\ cos(lus) + D\ cos(3a;s) + D l _ cos[(w - 2uj 8 )s] + D\ cos[(cj + 2cj e )s]. (61) 

The functions cos (lo s) and cos (tog s) represent a new problem for the third-order deviation, 
as they are resonance terms whose angular frequency lo (or ooq) is the same as the eigenvalue of 
the matrix-operator acting on the left-hand side, yielding secular terms, proportional to s. To 
avoid unbounded deviations, we can apply the Poincare's method |11] to take into account possible 
perturbation of the basic frequency itself, replacing lo (or loq) by an infinite series in powers of the 
infinitesimal parameter, which in our case can be the "eccentricity" e = 



LO 



UJr, 



00 



e ooi + e 2 co 2 + e 3 ^3 + • 



where the new lo is renamed oo p and loq is the old lo, and 

WO — ► 00Q p 



lo 90 + e toei + e 2 1^2 + e 3 Lo e3 + 



(62) 



(63) 



where the new tog is renamed loq v and o^g is the old loq. 

We shall build the complete differential equation for x^ 1 , taking together the harmonic oscillator 
equations for 5r, 5 2 r and 5 3 r: 



d 2 , r 8 2 r 5 3 r. 9/[ . 5 2 r 8 3 r. 



A 3 r A 3 ri . , A 2 r 2 , , 
— 1 — cos(u;ps) H — cos(2w p s) 
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A 3 r3 A 3 r_ A 3 r+ 
H — cos(3w p s) H - — cos[(w p - 2o; ep )s] H - — cos[(u;p + 2uq p )s\, (64) 

and for 59, 5 2 9 and 5 3 9: 

d 2 fsn 5 2 9 5 3 9^ 2 6 2 6 5 3 9^ A 2 6_ . . A 2 9 + 

j^K° e + — + -g~) + "oai™ + — + -g-) = cos[(ujp - uj0 P )s\ + — ^— cos[(w p + w e p)sj 

A 3 A 3 A 3 $ A 3 

H — cos(uJe P s) H — cos(3w ep s) H ^ cos[(2w p - w 0p )s] H cos[(2w p + oj0 p )s]. (65) 

Then, developing both sides into a series of powers of the parameter e, we can not only recover the 
former differential equations for the vectors 5 2 x fl , 5 3 x^, but get also some algebraic relations 
defining the corrections u\, u>2, ^ei an d ^92-, see the Appendix 2. In the Schwarzschild limit, we 
have: 

3 M 3/2 (6-37f) 

wi = 0, u 2 = -- nt - /9 , = • (66) 

4 i?5/2 /TTW (1 _ 6M )3 /2 



R R 

so the new frequency corrected by the Poincare's method is simply: 



_VM\ll~™ 3K) 2 M^ (6-37%) 



~ P i? 3 / 2 yJl-Uf 4 ^ 9/2 3M(1_ 6M)3/2' ^ 

Finally, we can obtain that the first and second-order deviations are the same, but with the new 
uj p and oj0 p ; and the third-order deviation 5 3 x^ is given by: 

5 3 r = 5 3 ro + 5 3 rs cos(3w p s) + <5 3 r_ cos[(cjp — 2uie p )s] + 5 3 r + cos[(w p + 2ujg p )s\, (68) 

5 3 9 = 5 3 9 3 cos(3ujg p s) + <5 3 6L cos[(2w p - uj 8p )s] + <5 3 (9 + cos[(2w p + uj 9p )s], (69) 

(5 3 </> = (5 3 4>o)s + 5 3 4>i sin(o;ps) + 5 3 <f>3 sin(3a;ps) + <5 3 </>_ sin[(w p — 2ujg p )s\ 

+ <5 3 + sin[(aj p + 2^ p )s], (70) 

S 3 t = (5 3 t )s + 5 3 h sin(o;ps) + 5 3 t 3 sin(3u;ps) + <5 3 i„ sm[(u p - 2uj ep )s] 

+ 5 3 t + sm[{u p + 2u 6p )s]. (71) 

We can choose the initial conditions of the differential solutions so that the constants 5 3 tq, 5 3 (f>o and 
5 3 to are null. The Appendix 2 lists the explicit values of the above coefficients. 

The long expressions of the Kerr case are well simplified in the Schwarzschild limit: 

S 3 9 = 0, (5 3 r_ = 5 3 r + = 5 3 <p„ = 5 3 (f> + = 5 3 t„ = 5 3 t + = 0, (72) 

and the non-null coefficients are: 

s s r _ -9K) 3 (2-28f +97^) 

3 8R 2 (1-^f) 2 ' [ ' 

3 , = 9(n 5 ) 3 K) 3 (26 ~ 336f + 1Q83^) 

91 i?3 (1 _6M)3/2' 93 4R 3 (1 _ 6M)5/2 ' ^ 

_ 3Kf [M (2- 19^+40^-36^) 

1_ i? 2 V« (i - ^) 3 (i - tt) 3/2 ' 
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3 4i?2 Y R (1 _ 2M )3(1 _ 6M)5/2 " 

The same approach could be used in the second-order deviation calculations, but it is not necessary 
because there are no resonances in the second-order deviation equations. 

Now we can compare the perihelion advance with the post-Newtonian limit, Eq. (jl|). The 
Schwarzschild limit gives the perihelion advance as 

[QttM 27itM 2 135ttM 3 \ (n r ) 2 f 9ttM 159ttM 2 

585 vrM 3 \ K) 4 f 81 vrM 2 594 vrM 3 \ 
+^3- + -J + S§h + + -J + - W 

which agrees with the Aip of Eq. (jlj) after replacing the major semi-axis a by the value of Eq. (|56|): 

6ttM 9vre 2 M . . 

A^ = — + — — + ... (78) 

The perihelion advance of Ay? in the Kerr case depends on M , R, a, Uq and n®, and can be described 
with high accuracy by a long explicit expression, not shown here. 

7 Discussion 

We have further developed a new method for calculating geodesies in a completely relativistic setting, 
using higher-order deviations without introducing the Newtonian or post-Newtonian approximations. 

The computation of first, second and third-order deviations for the Kerr metric have shown that 
this method can be reduced to a straightforward iteration of solving linear systems of differential 
equations with constant coefficients. 

The only complexity resides in the simplification of symbolic coefficients of the deviations, which 



is successfully performed by means of symbolic computing softwares [19]-[£(|. 

It is interesting to observe how at the very first level of approximation the angular momentum 
density a of the central body influences the perihelion advance via two different effects, which are 
linear and quadratic in a, respectively. The expressions linear in a depend on the sign of this 
parameter, i.e. on the relative sign of two rotations: that of the central body, and the direction of 
the orbital motion - a kind of spin-orbital coupling. This is the so-called dragging effect characteristic 
for General Relativity, which tends to raise the perihelion advance if the rotation of the planet is 
in the same direction as the rotation of the central body itself, and tends to decrease the perihelion 
advance if these two rotations are opposite to each other. The terms quadratic in a represent an 
additional perihelion advance which is due to the fact that the non-vanishing angular momentum of 
the central body is perceived from the exterior as an extra energy, which by the equivalence principle, 
may be considered as an extra mass 5M added to the central mass M; therefore, it always tends to 
produce higher perihelion advance. 

It is worth stressing that these effects are absent in the first-order post-Newtonian approximation. 
In this sense, our method gives a shorter way enabling one to display certain effects, than the 
commonly used post-Newtonian approach. Its convergence properties are very good, too, so that 
there are physical situations when it is more appropriate. Consider a small mass rotating quite 
close to a black hole, so that the quantity GM/rc 2 ~ v 2 /c 2 is of the order of 0.1; then the second 
post-Newtonian effects are of the order 0.01. Now, if the eccentricity of the orbit is of the same 
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order, i.e. e = 0.1 then our third-order terms give the precision of 0.001, keeping the quasi exact 
functional dependence on physical parameters GM/r and a. 

There are many possible applications and further developments. The computation of fourth and 
higher-order deviations in Schwarzschild and Kerr metrics can improve the accuracy for practical 
calculations, and is just a matter of spending more time and computer resources because we have 
developed a semi-automatic program for explicit calculation of higher-order geodesic deviations. 

The gravitomagnetic clock effect, i.e., the time-difference between the orbits of two freely counter- 
revolving test particles around a central rotating mass M, is an ideal target for high-order geodesic 



deviations, as the usual approaches are limited to circular orbits [21| or slowly rotating mass M [B^], 
i.e., small values of a. So the high-order geodesic deviations method has the potential to compute 
the gravitomagnetic clock effect in the case of strong general relativistic effects (large values of M 
and a). 

For practical applications involving the X-ray or gravitational radiation, the dynamics of accre- 
tion disks, etc, it would be useful to generalize our method to the case of initial orbits inclined w.r.t. 
the equatorial plane of the rotating central mass M. For example, the Ref. p^J considers inclined 
orbits in the Kerr metric, with the constraint of low-eccentricity orbits, i.e., up to the first-order 
deviation. 

Still within the test particle concept, we can extend it for test bodies carrying charge and/or in- 
ternal spin, see Ref. [24] where the first-order geodesic deviation is derived in the Reissner-Nordstr0m 



background field. We foresee that higher-order geodesic deviations for test particles with spin can 
provide useful results to compare with the experimental data of satellite gyroscopes. 

We can also replace the background metric by some cylindrical or axially-symmetric metric to 
investigate approximated models of star and galaxy orbits, accretion disks, etc. 

Almost all the work available in the domain of gravitational radiation is based on post-Newtonian 
approximations |25|-p9|1. Within the higher-order geodesic deviations approach, one possibility 
is to maintain the test particle mass m negligible compared to M, and compute the emission of 



gravitational radiation [30] with some formula better suited for Schwarzschild and Kerr metrics than 
the quadrupole formula [|31j] . Another possibility, more challenging, is to cope with the finite-size of 
the mass m by taking it into account with appropriate perturbation of the background metric, then 
trying to repeat the higher-order geodesic deviations calculations and finally employing a modified 
gravitational radiation formula based on the perturbed metric. 
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Appendix 1 

The coefficients of the solution, in the Kerr metric case, for the second-order geodesic deviation 
5 2 x ll (s) are: 

Pe- = -Pe+ = ^^^, (so) 

R v h 
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Appendix 2 

With the third-order deviation in the Kerr metric, we can determine the values of: 

LOl = LO01 = 0, 
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therefore the new frequencies corrected by the Poincare's method, which are exact up to the second 
order w.r.t. the small parameter e = are given by: 



w + 



^2 , ^>0p = ^00 + 



^02 



(86) 



The coefficients of the solution, in the Kerr metric case, for the third-order geodesic deviation 
S 3 x^(s) are: 
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